#clean up
rm(list=ls())

#options
options(scipen=10)

#load libraries
library(foreign)
library(MCMCpack)
library(lmtest)
library(car)
library(xtable)

#set working directory
#setwd('/Users/jamie/Dropbox/hedged lobbying/')
#setwd('/Volumes/MONOGAN/Hedged Lobbying/')

#load data
hedge<-read.dta('lobbyHedgeRep.dta')

#Original OLS model
ols.orig<-lm(Dtotallobbymoney~DAUMHedgeFunds+DGDP+Dbudget,data=hedge);summary(ols.orig)
bgtest(ols.orig) 
bptest(ols.orig,studentize=FALSE)
vif(ols.orig)
influenceIndexPlot(ols.orig, vars=c("Cook","Studentized","hat"), id.n=2, labels=hedge$year[-1]) #2008 is huge on Cook's D
plot(y=ols.orig$residuals,x=ols.orig$fitted.values)
plot(y=ols.orig$residuals,x=hedge$year[-1])
plot(y=ols.orig$residuals,x=hedge$DAUMHedgeFunds[-1])
plot(y=ols.orig$residuals,x=hedge$DGDP[-1])

#Correcting for 2008 being an influential data point
ols.fin<-lm(Dtotallobbymoney~DAUMHedgeFunds+DGDP+Dbudget+var2008,data=hedge);summary(ols.fin)
bgtest(ols.fin)
bptest(ols.fin,studentize=FALSE)
vif(ols.fin)
influenceIndexPlot(ols.fin, vars=c("Cook","Studentized","hat"), id.n=2, labels=hedge$year[-1])#no wild Cook's D values
plot(y=ols.fin$residuals,x=ols.fin$fitted.values)
plot(y=ols.fin$residuals,x=hedge$year[-1])
plot(y=ols.fin$residuals,x=hedge$DAUMHedgeFunds[-1])
plot(y=ols.fin$residuals,x=hedge$DGDP[-1])

#Bayesian version of the final model, 2008 as an influence point
mcmc.fin<-MCMCregress(Dtotallobbymoney~DAUMHedgeFunds+DGDP+Dbudget+var2008,data=hedge,burnin=10000,mcmc=100000,b0=0,B0=.001,marginal.likelihood="Chib95")
summary(mcmc.fin, quantiles=c(.05,.10,.5,.90,.95))
geweke.diag(mcmc.fin, frac1=0.1, frac2=0.5)

###WRITE OUTPUT FOR FINAL MODEL###
#TABLE
output<-summary(mcmc.fin, quantiles=c(.05,.10,.5,.90,.95))
output.x<-xtable(cbind(output$statistics[,1:2],output$quantiles[,c(1,5)]),digits=4)
print(output.x)

#FIGURES
#postscript('hedgeSeries.eps', paper='special', horizontal=FALSE, onefile=FALSE, width=6, height=6, pointsize=10)
#png('hedgeSeries.png',width=1800,height=1800,pointsize=31)
par(mfcol=c(2,2))
plot(y=hedge$totallobbymoney,x=hedge$year,ylab="Total Money on Lobbying (Billions)",xlab="Year",type="l",main="Lobbying Spending")
plot(y=hedge$GDP,x=hedge$year,ylab="GDP (Trillions)",xlab="Year",type="l",main="U.S. GDP")
plot(y=hedge$AUMHedgeFunds,x=hedge$year,ylab="Hedge Fund Assets (Billions)",xlab="Year",type="l",main="Hedge Fund AUM")
plot(y=hedge$budget,x=hedge$year,ylab="Federal Outlays (Trillions)",xlab="Year",type="l",main="Federal Budget")
#dev.off()

#postscript('hedgeDensities.eps', paper='special', horizontal=FALSE, onefile=FALSE, width=6, height=3, pointsize=10)
#png('hedgeDensities.png',width=1800,height=900,pointsize=30)
par(mfcol=c(1,2))
quantile(mcmc.fin[,2],.039) #positive, 96.1% probability
dens.2<-density(mcmc.fin[,2])
plot(dens.2,xlab="Coefficient for Change in Assets Under Management",ylab="Density",main="")
qLOW <- quantile(mcmc.fin[,2],.039)
q100 <- quantile(mcmc.fin[,2], 1)
x1<-min(which(dens.2$x>=qLOW))
x2<-max(which(dens.2$x<q100))
with(dens.2, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col="gray80",border=NA))
text(x=.0005,y=250,"p=0.961")

quantile(mcmc.fin[,4],.034) #positive, 96.6% probability
dens.3<-density(mcmc.fin[,4])
plot(dens.3,xlab="Coefficient for Change in Federal Budget",ylab="Density",main="")
qLOW <- quantile(mcmc.fin[,4],.034)
q100 <- quantile(mcmc.fin[,4], 1)
x1<-min(which(dens.3$x>=qLOW))
x2<-max(which(dens.3$x<q100))
with(dens.3, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col="gray80",border=NA))
text(x=.5,y=.25,"p=0.966")
#dev.off()


